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Abstract 

In this work a method is proposed for computing step free energies for faceted solid-liquid interfaces 
based on atomistic simulations. The method is demonstrated in an application to (111) interfaces in ele- 
mental Si, modeled with the classical Stillinger- Weber potential. The approach makes use of an adiabatic 
trapping procedure, and involves simulations of systems with coexisting solid and liquid phases separated 
by faceted interfaces containing islands with different sizes, for which the corresponding equilibrium 
temperatures are computed. We demonstrate that the calculated coexistence temperature is strongly 
affected by the geometry of the interface. We find that island radius is inversely proportional to super- 
heating, allowing us to compute the step free energy by fitting simulation data within the formalism of 
classical nucleation theory. The step free energy value is computed to be j st = 0.103 ± 0.005 x 10~ 10 
J/m. The approach outlined in this work paves the way to the calculation of step free energies relevant 
to the solidification of faceted crystals from liquid mixtures, as encountered in nanowire growth by the 
vapor-liquid-solid mechanism and in alloy casting. The present work also shows that at low undercool- 
ings the Stillinger- Weber interatomic potential for Si tends to crystallize in the wurtzite, rather than the 
diamond-cubic structure. 
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I. INTRODUCTION 



The properties of solid-liquid interfaces are known to play critical roles in governing defect 
formation and the evolution of micro structural patterns during the solidification of a crystal 
from its melt. 1 In applications of crystal- growth theories and mesoscale simulation methods to 
the study of solidification phenomena, quantitative information is often required for the ther- 
modynamic and kinetic properties of solid-liquid interfaces, and the degrees to which these 
properties vary with crystal orientation. For materials such as metals where crystal-melt inter- 
faces are molecularly rough, the most relevant properties in this context include the magnitude 
and crystalline anisotropy of the interfacial stiffness and mobility,^ which dictate, respectively, 
the capillary undercooling at a growing crystal-melt interface, and the relationship between lo- 
cal undercooling and interface velocity. For many crystalline compounds and covalent solids, 
the degree of crystalline anisotropy in crystal-melt interfacial free energies can be highly pro- 
nounced, such that the equilibrium crystal Wulff shape is characterized by faceted orientations. 
In this case, crystal-growth kinetics are typically strongly influenced by the properties of steps 
at faceted solid-liquid interfacesP^For example, the excess free energies and mobilities of these 
linear defects play critical roles in governing the rate of nucleation and the growth of crystal is- 
lands during layer-by-layer growth of faceted crystals, as observed in the synthesis of nanowires 
by the vapor-liquid-solid mechanismP^ 

Due to the difficulty inherent in performing direct experimental measurements of the 
thermodynamic and kinetic properties of crystal-melt interfaces, numerous computational 
approaches have been developed for their calculation within the framework of atomistic 
simulations.^ Specifically, for crystal-melt interfaces in systems with molecularly-rough inter- 
faces, equilibriurrP^G^and non-equilibriurrP^molecular-dynamics (MD) and metadynamics^^ 
methods have been developed for the calculation of interfacial free energies and associated 
crystalline anisotropics. Similarly, both equilibrium and non-equilibrium MD methods have 
been applied in calculations of the magnitude and anisotropics of the mobilities of such 
interfaces] 1 * 25 ^ 26 ! To date, far less simulation work has been devoted to the investigation of 
faceted crystal-melt interfaces.^^ Additionally, the properties of steps at faceted crystal-melt 
interfaces have been investigated primarily by kinetic Monte-Carlo simulations 33 to date. To the 
best of the authors' knowledge, only one MD study has been applied to the calculation of step 
mobilities, 29 and no direct MD calculations of step free energies have yet been published. As a 
consequence, the molecular-level understanding of the properties of steps at faceted solid-liquid 
interfaces remains less advanced relative to those for rough crystal-melt interfaces or steps at 
crystalline surfaces. 
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In this work we propose a method for the calculation of step free energies 7 si at faceted 
crystal-melt interfaces from equilibrium MD simulations. The approach makes use of simula- 
tion geometries involving coexisting solid and liquid phases, separated by a faceted solid-liquid 
interface containing an island of solid or liquid on an otherwise flat terrace. As described in 
the next section, the coexistence temperature measured in such simulations is expected to vary 
with the curvature (radius) of the islands, such that the relationship between radius and the co- 
existence temperature measured in MD can be used to extract orientation- averaged values of 
7 S (. The approach is illustrated in this work in the application to steps in the Stillinger- Weber 
model of SiP^The simulation methodology employed in this application is described in Section 



III, followed by a presentation of results and a summary and discussion in Sections IV andfV] 



respectively. In Section [IV] we also present results illustrating that the Stillinger- Weber model 
of Si shows a tendency to crystallize in the wurtzite, rather than the diamond-cubic, structure at 
relatively low undercoolings. While this finding is not directly relevant to the present study, it is 
important for the detailed interpretation of simulations of crystallization based on this potential. 



II. THEORETICAL FORMALISM 

In Fig. [T] we consider three possible geometries of a faceted solid-liquid interface, where 
terraces of two different "heights" are present, separated by steps. In Figs [T]a and c the ge- 
ometry corresponds to circular island of solid and liquid phases, respectively. Figure [T]b by 
contrast shows a situation where terraces of different height are separated by straight steps. 
When amounts of solid and liquid phases change due to melting or crystallization the length of 
the step boundary changes in systems with circular geometry, while it remains constant for the 
planar case. Thus, equilibrium conditions for planar and circular geometries should be different. 

We consider in more detail the equilibrium conditions for the systems with a circular island, 
using the framework of classical nucleation theory (CNT).^Sl Consider a circular solid island 
on a liquid terrace. The free energy of the system relative to a state without an island is given 
by 

G = -A[ip A irr 2 + 2Trrj st , (1) 

where A/i = hi — n s is the difference between the chemical potential in the liquid and solid 
phases, pa is the atomic density per unit area of the solid plane that forms the island, and r is 
the radius of the island. The size of the critical nucleus R corresponds to a maximum of G(r) 
and is given by: 
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jj> _ ^stTm ^2~j 

pAH m {T m — T) 

where we have used Afj, = ^P- (T m — T), with T m corresponding to the equilibrium temperature 
for a planar step geometry and H m is a heat of melting per atompEO According to Eq. §2§ a 
system with a solid island of radius R is in equilibrium at temperature T given as: 

T = T m (l - d/R) (3) 
where the capillary length d is given as: 

d = (4) 

Note that, from Eqs. and Q a solid island has a coexistence temperature below the 
bulk melting point. By contrast a liquid island has a coexistence temperature above T m and 
corresponds to a superheated state. Finally, for the interface geometry with a planar step shown 
in Fig [TJb the coexistence temperature is simply T m , which follows from Eq. ([3]) by taking 
the limit R — >■ oo. Clearly, in the latter case the equilibrium temperature is independent of 
island size. The formalism discussed above neglects curvature corrections to 7 st as well as 
step interactions and stress effects. However, these effects can in principle be incorporated by 
appropriate generalizations of the thermodynamic formalism. 

The MD approach considered in this work makes use of the thermodynamic formalism de- 
scribed in this section as follows. By performing equilibrium MD simulations in the micro- 
canonical ensemble (NVE), the constraint of fixed total system energy can be used to stabilize 
the various types of equilibrium states described above. By changing the system energy, islands 
of different radii can be stabilized and the corresponding system temperature derived. From the 
relationship between R and T, the step free energy can then be determined. 



III. METHODOLOGY OF ATOMISTIC SIMULATIONS 

We demonstrate the approach outlined in the previous section by considering the Stillinger- 
Weber model of elemental SiP^ This model represents one of the simplest classical inter- 
atomic potentials that gives rise to faceted solid-liquid interfaces ) 27 l 37 l 38 l Several different val- 
ues of the melting temperature have been reported for this potential, from simulations with 
coexisting solids and liquids separated by (111) oriented interfaces, namely: T m = 1677-fT,^ 
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T m = 1685K 38 and T m = 1 69 IK*® Unfortunately, it was not possible from these studies to 
know whether the coexistence conditions corresponded to any of the different types of interface 
geometries described in Fig. [T] or to a flat terrace with no steps. We note that, in the latter case, 
the coexistence temperature derived in an MD simulation can be different from the true melting 
point, since the system can spend a significant amount of time in such a state, 4 " due to the free 
energy barrier required to nucleate islands of solid or liquid. During this time thermal (uniform 
temperature) and mechanical equilibrium between the phases can be established, but the phases 
are not in true equilibrium with respect to phase change fluctuations, i.e. the chemical potentials 
of the solid and liquid phases are not equal. To enable phase fluctuations on the time scale of 
MD simulations it is desirable to have preexisting steps at the interface. All equilibrium states 
implemented in this study satisfy this requirement. 

A. Simulation block 

To apply the simulation approach described above, we began by creating a simulation block 
with diamond-cubic-structured solid and liquid phases as shown schematically in Fig. [2}a. The 
(lll)-oriented solid-liquid interface was perpendicular to the z direction of the block. The x 
and y directions were parallel to crystallographic directions [110] and [112] of the solid phase, 
respectively. The phases were equilibrated for several nanoseconds in a microcanonical (NVE) 
ensemble. In the resulting equilibrated state the dimensions of each phase were approximately 
cubic. In the next step, we selected a region containing the solid-liquid interface as well as 
regions of homogeneous solid and liquid phases as shown in Fig. [2}b. Atoms inside the top and 
bottom layers were designated as boundary regions. The thickness of each boundary layer was 
10 A. 

This procedure was used to create several simulation blocks with different sizes. The number 
of atoms in the resulting systems ranged from 25196 to 226744. The dimensions L x and L y 
parallel to the interface ranged from 5 to 40 nm. The dimension normal to the interface L z =4.5 
nm was the same in all blocks. Periodic boundary conditions were applied in the x and y 
directions. In the z directions the boundary conditions were not periodic. The positions of 
atoms inside the solid boundary region were fixed during the simulations. The top boundary 
region of liquid atoms moved as a rigid body during the simulations to ensure zero pressure 
inside the liquid phase below the boundary region. 

Such boundary conditions were chosen for two reasons. First, to apply the formalism de- 
scribed in the previous section it is necessary to model only one solid-liquid interface, which is 
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impossible with periodic boundary conditions in the z direction. Second, the rigid body bound- 
ary condition was chosen over the open liquid surface to isolate the solid-liquid interface from 
the effect of surface capillary waves and density variations. Thus, the rigid slab of liquid (and 
solid) mimics an infinitely large homogeneous phase. Due to the use of relatively small lengths 
in z the values of the step energies computed below could be affected by finite-size effects, 
because we are constraining the phonon spectrum in the crystal and the density oscillations in 
the liquid near the solid-liquid interface. We emphasize however that the primary purpose of 
the present work is to demonstrate the adiabatic-trapping methodology for computing step free 
energies. Finite size effects could be straightforwardly investigated in future studies focused on 
precise values for specific systems. 

B. MD simulations 

Molecular dynamics simulations were performed using the Large-scale Atomic/Molecular 
Massively Parallel Simulator (LAMMPS) software package.^ Thermal expansion of solid Si at 
T = 1679 K was computed from a simulation in an NPT ensemble with zero pressure. Constant 
temperature and pressure in the simulation were imposed using a Nose-Hoover thermostat an d 
an Anderson barostat, 43 respectively. 

Molecular Dynamic simulations of the solid-liquid interface were performed in a micro- 
canonical ensemble (NVE) for times up to 50 ns after 5ns equilibration stage, using a time step 
of 1 fs. During the simulations x and y dimensions of the simulation block were kept fixed, 
with the lattice parameter in the solid part of the block equal to the stress free lattice parameter 
at T = 1679 K. 

During the simulations snapshots of the system were saved every 0.5 ns. The snapshots 
contained positions of atoms, energies and stresses and were used for post processing. The 
equilibrium temperature T of the simulation was calculated from an average of the kinetic 
energy over the production stage. 

1. Trapping procedure 

In canonical simulations (NVT) solid/liquid islands on liquid/solid terraces are usually not in 
equilibrium with the rest of the system: subcritical islands shrink, whereas supercritical islands 
grow to completely cover the terrace. Critical islands correspond to a state of unstable equilib- 
rium. To compute sizes of critical nuclei we employed a trapping procedure used previously 
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in a study of heterogeneous nucleation at grain boundaries. 17 In a finite system and a restricted 
ensemble (NVE) equilibrium between the island and the rest of the system can be stabilized due 
to the constraint of constant energy. To achieve a state of a stable equilibrium, fluctuations in 
the size of the island must produce sufficient changes in temperature of the system. As a con- 
sequence the z dimension of the simulation block cannot be too large. Also lateral dimensions 
of the island have to be comparable to x and y dimensions of the simulation block. For the 
reasons above, for each size of the simulation block a structure with an island of a size within 
a certain range had to be carefully prepared. If the size of an island was too small, it could 
completely disappear by a fluctuation during a 50 ns long simulation. Additionally, a very large 
island could hit the boundary conditions and change its morphology from circular to planar, 
thus changing the nature of the equilibrium state. 

To prepare a simulation block with a suitable size of an island, the radius of the solid/liquid 
island on the terrace can be modified by injecting/removing heat into a system. This was 
achieved by modifying velocities of atoms in the input file and equilibrating the system during a 
several nanosecond run in the NVE ensemble. We prepared simulation blocks with solid/liquid 
islands embedded inside the liquid/solid terrace as well as blocks with solid and liquid terraces 
separated by a planar step boundary. The former states correspond to curved equilibrium with 
equilibrium temperature depending on the size R of the islands, while the later case corresponds 
to planar equilibrium with temperature independent of the sizes of the terraces. 

2. Data analysis 

To analyze the structure of the interface and calculate the sizes of islands it was necessary 
to identify atoms belonging to solid or liquid phases. This was done using a structural order 
parameter employed in the previous studies of the Stillinger- Weber Si systemP^H The calcu- 
lated order parameter takes values close to 1 inside solid phase and inside liquid phase. A 
threshold value of the order parameter was calculated as an average of the values inside solid 
and liquid phases. Atoms with the value of the order parameter higher than this threshold value 
were identified as solid, while all other atoms were identified as liquid. 

The structure of the solid-liquid interface was analyzed by identifying isolated solid/liquid 
islands on the liquid/solid terrace. For each atom belonging to solid/liquid phase a list of the 
nearest neighboring atoms of the same phase were constructed. Then a random walk through 
neighboring atoms was implemented in the following way. First, a random atom of a given 
phase was selected and assigned a cluster number. Second, its neighbor belonging to the same 
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phase was randomly selected and assigned the same cluster number. Thus, repeating the pro- 
cedure all atoms belonging to one cluster can be identified. The procedure was previously 
employed to identify clusters in a binary liquid.^ The largest island identified represents the 
nucleus, while all other islands are fluctuations containing only a few atoms. These small is- 
lands are present in the structure of the terrace at a given temperature.^ 

Our cluster analysis of the interface allows us to calculate the center of mass of the largest 
island on the terrace. To compute an average size and shape of the island we average the atomic 
density of the instantaneous islands relative to instantaneous centers of mass. Once the average 
atomic density is computed, the radius R of an island is calculated as a distance from the center 
of mass to the point where the average density value is equal half of the density of solid (1 11) 
layer. 

IV. RESULTS 

A. Interface Geometries 

Fig. [3] shows a cross-sectional view of a region of a (111) solid-liquid interface containing 
a step separating two terraces of different heights. Atoms belonging to the liquid phase are 
colored in red, while solid atoms appear in blue. The [110] crystallographic direction in the 
solid phase is normal to the plane of the figure. The top solid (111) plane is incomplete and 
forms a step boundary with the liquid terrace. The height of the step is equal to the bilayer 
spacing of the diamond-cubic structure along the [111] direction. 

Fig. [4] shows plan views of snapshots from the equilibrium configurations of faceted (111) 
solid-liquid interfaces with different island geometries. Figs[4ja and b represent curved equilib- 
rium states of a solid island on a liquid terrace, and liquid island on a solid terrace, respectively. 
Both simulation blocks have the same size with lateral dimensions of 30 x 30 nm. These snap- 
shots illustrate that the instantaneous structure of the step is extremely rough, with the shape of 
the island being quite complicated and characterized by "overhangs" and large deviations from 
a perfect circular geometry. The equilibrium temperatures calculated for these solid and liquid 
islands were T ~ 1678 K and T = 1685.4 ± 0.7 K, respectively. Thus, for the same size of the 
simulation block the calculated equilibrium temperature was not the same due to the difference 
in geometry of the atomic layer at the solid-liquid interface. 

Fig. |4]c shows a plan-view snapshot illustrating the equilibrium state of the interface with 
a planar geometry of the steps at the calculated temperature T = 1682.4 ± 0.6 K. The z di- 
mension normal to the solid-liquid interface was the same as in all other blocks studied. The 
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lateral dimensions were 10 x 40 nm 2 parallel to the x and y axis respectively. The simulation 
block contains two step boundaries. This geometry with one dimension of the simulation block 
smaller than the other is more advantages for simulations of planar-step equilibrium. It ensures 
that the liquid or solid part of the interface layer does not transform into an isolated island with 
circular boundaries during a long simulation. Simulations with planar step geometries were also 
performed for systems with lateral periodic dimensions of 5 x 10 nm 2 and 10 x 30 nm 2 . Identi- 
cal equilibrium temperatures were calculated for these blocks within the statistical accuracy of 
the simulations, as discussed in further detail below. 

1. Wurtzite Formation 

In the simulations with solid circular islands in a liquid terrace we observed nucleation of 
new islands showing a registry with the underlying solid layer corresponding to the wurtzite, 
rather than the diamond cubic structure. These new "wurtzite" islands were observed to grow 
and consume the preexisting "diamond" islands. After this transformation the islands were 
observed to remain in the wurtzite state. Previous MD simulations of solidification with (111) 
oriented solid-liquid interfaces 38 and of isolated solid nuclei 45 reported that the Stillinger- Weber 
Si potential solidifies in a random mixture of stacking sequences. These simulations modeled 
large undercoolings, and the observations could be a consequence of rapid solidification ki- 
netics. To investigate this issue further, we performed a separate simulation of solidification 
at 1670 K (-12K undercooling) in the NP 2 T ensemble, featuring fully periodic boundary con- 
ditions with fixed area parallel to the interface and dynamic periodic lengths normal to the 
interface to ensure zero normal stress. The simulation block had dimensions 5 x 5 x 18 nm 3 , 
with two (1 1 1) solid-liquid interfaces normal to the z direction. Figs. [5Ja and b show the initial 
state of the simulation block with the solid phase in the diamond structure and the simulation 
block after 110 ns long simulation, respectively. The liquid phase crystallized into wurtzite 
phase. Only one out of 32 newly grown (111) layers was observed to have the diamond-cubic 
stacking sequence. 

The observation of the nucleation and growth of the wurtzite structure during simulations 
at low undercoolings is likely an artifact of the Stillinger- Weber potential. Specifically, due to 
its short cutoff distance the potential gives identical zero-temperature energies for the diamond 
and wurtzite structures. 45 The wurtzite phase may appear during solidification for two reasons. 
First, it may be more stable than the diamond phase at high temperatures, which means that 
it would have a higher melting point. Second, the liquid may crystallize into wurtzite due to 
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purely kinetic effects. For example, if the step free energy of a wurtzite island is lower than 
that of a diamond island, a monolayer of wurtzite will nucleate first. At high undercoolings 
nucleation barriers for diamond and wurtzite structures become nearly zero and the nucleated 
structure is determined by random fluctuations. As a result, it is understandable that at large 
undercoolings random sequence of these two phases can be observed. 38 The melting point of 
the wurtzite phase was calculated using simulation blocks with planar geometries for the steps 
and was found to be TjJ = 1681.5 ± 1.3 K, which is nearly identical with that of the diamond 
phase. Therefore, the second scenario, involving kinetic stabilization of the wurtzite phase, is 
more consistent with our simulation results. 

B. Island Radius versus Temperature and Step Free Energy 

Because of the diamond to wurtzite transformation we were unable to calculate the rela- 
tionship between island radius (R) and temperature (T) for diamond-cubic solid islands on the 
liquid terrace. Accurate equilibrium temperature calculations require tens of nano-seconds of 
sampling time and we were unable to accumulate the necessary statistics before transformations 
to the wurtzite structure occurred. In the simulations with liquid islands on solid terrace, nucle- 
ation of wurtzite islands is unfavorable because the temperature of the simulation is higher then 
the meting point of the wurtzite phase. As a result we were able to model liquid islands for any 
desired period of time to compute the relationship between R and T. 

Fig. [^demonstrates an average equilibrium shape of a liquid island on solid terrace at T = 
1685.4 ± 0.7 K. The shading in the image corresponds to the average density of solid atoms, 
as indicated in the legend. The solid-atom density inside the island is nearly zero, while it is 
close to the areal density of a (111) solid plane away from the island. The thickness of the 
transition region where the density changes between these limiting values is about 5 nm, which 
is comparable with the radius of the island in Fig. [6} While the instantaneous shapes of the 
island are very complicated (c.f., Fig. [4] a and b) the average shape shown in Fig. [6] appears to 
be almost perfectly circular. Roughness of the step as well as instantaneous deviations from the 
circular shape and size fluctuations all contribute to the broadening of this region. While the 
first two processes are most likely to be the dominant ones, the third process is always present 
and strongly affected by the size of the block and the island. 

Fig. [7] summarizes calculations of the equilibrium island radius and corresponding equilib- 
rium temperatures. The discrete points on the figure correspond to individual MD simulations. 
Squares correspond to simulations having steps with planar geometry (e.g., Fig. |4}c), while 
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other symbols represent simulations of isolated liquid islands on a solid terrace (e.g., the ge- 
ometry shown in Figs. [4|b and [6]). The continuous lines on the plot were obtained by fitting 
the curved equilibrium data points with Eq. ([2]) using 7 st and T m as fitting parameters. The 
data points with planar step geometry were not included in the fitting. The excellent agree- 
ment between the fit and the discrete data points allows us to conclude that equilibrium island 
radius R is inversely proportional to AT = T — T m . The melting point obtained from the fit 
was T m = 1681.96 ± 0.05 K, which is identical within statistical uncertainties with the equilib- 
rium melting temperature calculated directly from MD simulations with planar step geometries: 
^pLnar = 1682.0 ± 1.0 K, computed as an average over the three independent simulations with 
different sizes. 

Using the data for latent heat of melting for Stillinger- Weber Si, H m = 31 kJ/mole (5 x 10~ 20 
J/atom)j22 the areal atomic density of the (111) plane = 0.157 x 10 20 m~ 2 , and the value 
of the melting point given in the previous paragraph, we calculate from the fit shown in Fig. [7] 
a value for the step free energy of 7 st = 0.103 ± 0.005 x 10~ 10 J/m, where the error bar was 
determined from the fitting. Dividing this value by the (111) interplanar distance S 111 ^ = 3.14 
A, we obtain a value for the so-called perimeter free energy® of 7^ = 0.033±0.0016 J/m 2 . The 
perimeter free energy has units of energy per unit area and can be compared with solid-liquid 
interface free energy. 

The solid-liquid interface free energies for Stillinger- Weber Si calculated using the cleaving 
technique 33 for (100), (111) and (110) orientations are 0.42 ± 0.02 , 0.34 ± 0.02 and 0.35 ± 
0.03 J/m 2 respectively. These values are approximately an order of magnitude larger than the 
perimeter free energy quoted above. Similarly, multiplying the interface free energy by the 
thickness of one (111) monolayer would give an estimate for 7 si of 1.16 x 10 10 J/m, which is 
an order of magnitude larger than the step free energy calculated by MD in the present work. 

V. DISCUSSION AND CONCLUSIONS 

In this work we presented a method for MD calculations of the free energies of steps of 
faceted solid-liquid interfaces. The approach makes use of simulations in which the faceted 
interface contains islands with curved steps of differing radii. For a range of island radii the 
coexistence temperature is derived from equilibrium MD simulations and the resulting relation- 
ship between R and T — T m is used to back out 7^. 

The approach was demonstrated in an application to (1 1 1) faceted interfaces in the Stillinger- 
Weber model of elemental Si. We simulated interfaces with planar geometries of the steps 
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as well as shapes of islands that were circular on average. In the later case, there were two 
different states, corresponding to a solid island in a liquid terrace, or a liquid island in a solid 
terrace. We demonstrated that the calculated equilibrium temperature depends on structure 
of the interface. When the step separating solid and liquid parts of the interface is planar, the 
calculated equilibrium temperature T m = 1682. 0± 1.0 K was found to be remarkably consistent 
and independent of the size of the islands. In the case of equilibrium with a circular island 
geometry, the island radius R was found to be inversely proportional to T — T m . 

The island size dependence on equilibrium temperature allowed us to compute step free 
energy 7^. The calculated step free energy is an order of magnitude smaller than the product 
of interface free energy and thickness of a (111) monolayer. The difference between the two 
quantities is not surprising. When a faceted interface undergoes a roughening transition, the 
step free energy becomes zero, while solid-liquid interface free energy generally remains finite. 
Therefore, these two properties are not directly related. Large island shape fluctuations observed 
in MD and the roughness of the step structures correlate with small value of the step free energy, 
and indicates that the interface may not be far from a roughening transition. 

A relatively wide range of equilibrium temperatures have been reported in the literature 
for coexistence simulations with (111) oriented interface. Our simulations provide an insight 
why, in general, this should be the case. The equilibrium temperatures calculated in this work 
ranged from 1682 K to 1695 K, depending on the exact nature of the interface geometry. This 
range would actually be twice large, if we include the undercooled branch where quantitative 
analysis of the data was not possible due to the formation of layers with wurtzite stacking. The 
current results thus illustrate that variations of equilibrium temperature within 26 K range are 
possible when simulating solid-liquid systems separated by the faceted (111) interface. During 
solid-liquid coexistence simulations with (111) interfaces with large bulk phases, the interface 
geometry may transition by fluctuations between undercooled, planar or superheated states. The 
interface can also adopt a geometry when the interface layer is complete (state with no step) 
and spend a significant amount of time trapped in this metastable state. As a result, depending 
on the initial conditions, amounts of phases and time of calculations, the resulting average 
temperature may deviate significantly from the true melting temperature. This study shows 
that for a consistent calculation of T m using the coexistence approach, for systems with faceted 
interfaces, an interface with a planar geometry of a step is desirable. 

The approach of modeling of plane and curved equilibrium islands to extract values of step 
free energies can be extended to multicomponent systems. In the later case MD simulations 
in an NVT ensemble or constant composition Monte CarlcP' simulations are more appropri- 
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ate. The stability of the island is maintained by constant average composition in the system, 
rather than a constant energy constraint. Performing a series of isothermal calculations at dif- 
ferent temperatures, one could recover the dependence of the step free energy on composi- 
tion/temperature. This dependence is an important ingredient in understanding and mesoscale 
modeling of layer-by-layer growth in the context of important technological processes such as 
nanowire growth by the vapor-liquid-solid mechanism. For example, in applications of such an 
approach to Si in contact with a eutectic Au-Si melt, the equilibrium temperatures probed in the 
simulations would be significantly lower than melting point of elemental Si. It is possible, that 
the step structure will not be as rough as than observed in the present simulation. Specifically, 
it is possible that the steps would become significantly straighter, such that the islands would 
develop "polygonal" shapes, using the terminology from Ref. ®. Investigation of equilibrium 
shape of such islands and their evolution during growth would be of fundamental interest. 

Finally, this work points to an important shortcoming of the Stillinger- Weber potential in 
modeling crystal growth in Si. Due to the small radius of interaction, the potential does not dif- 
ferentiate between the diamond and the wurtzite structures, and the energies of these structures 
are identical at zero temperature. 45 In the present simulations we observed that at low under- 
coolings the liquid crystallizes into the wurtzite structure, and interfacial islands with diamond 
stacking are unstable with respect to growth of islands with wurtzite stacking. The appearance 
of the wurtzite phase is often observed during nano-wire growth 49 and can be due to a difference 
in step free energies. It is possible that at different undercoolings the liquid would crystallize 
into phases are that are metastable in the bulk for kinetic reasons. To explore this issue further 
for growth of Si crystals would be of great interest, but for such purposes it would be desirable 
to employ a potential that better describes the relative energetics of the competing diamond and 
wurtzite phases. 
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Figure 1. Possible geometries of the solid-liquid interface with steps, a) Circular solid island on liquid 
terrace, b) Equilibrium with planar geometry of the step corresponds to the melting point T m . c) Circular 
liquid island on solid terrace. 
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Figure 2. Construction of the simulation block for modeling of an isolated interface, a) Bulk solid and 
liquid at equilibrium. Dimensions of each phase are approximately cubic, b) A region near interface is 
selected with dimension normal the interface significantly smaller than the lateral dimensions to create a 
new simulation block, c) A new simulation block with two boundary regions. The solid boundary region 
is fixed in subsequent simulations. The liquid boundary region moves as a rigid body to ensure zero 
pressure in the liquid. 
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Figure 3. (Ill) solid-liquid interface with a step. Solid atoms are colored in blue while liquid atoms are 
shown in red. The coloring is according to the order parameter described in the main text. The upper 
(111) plane is incomplete and forms a step with the liquid terrace on the right. The image was produced 
with the ATOMEYE visualization programpS 
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Figure 4. Structures of the solid liquid-interfaces with different types of terrace geometries, (a) A 
snapshot of a solid island inside a liquid terrace at T ~ 1678 K (undercooling), (b) Snapshot of a liquid 
island inside a solid terrace at T = 1685.4 ± 0.7 K (superheating), (a) and (b) correspond to curved 
equilibrium, (c) Solid and liquid terraces separated by a planar step boundary at T m = 1682.0 i 1.0 K. 
The images were produced with the ATOMEYE visualization program.^ 
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Figure 5. Simulations of solidification at T = 1670 K. Initially the simulation block contains liquid 
and a crystalline solid in the diamond structure. The liquid phase in a) solidifies primarily into wurtzite 
structure in b). The images were produced with the ATOMEYE visualization program.^ 
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Figure 6. Average density of atoms identified as solid according to the order parameter at T m = 1685. 4± 
0.7 K obtained by averaging of instantaneous snapshots. The orange region corresponds to a solid terrace 
with density close to that of a perfect (111) plane (pa = 0.157 x 10~ 20 m~ 2 ). The near zero density 
of solid atoms in the dark middle region of the plot represents the liquid nucleus. This density map was 
used to compute the size of the critical nucleus R. 
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Figure 7. Radius of critical nucleus ( o, A ) as a function of temperature. The dimensions of the 
simulation blocks parallel to the solid-liquid interface are indicated on the plot. The melting point for 
planar step geometries was calculated for different sizes of the simulation block ( □ ) . The continuous 
line was fit to discrete data points using Eq. Q. The dashed line indicates the melting point obtained 
from fitting the curved equilibrium points with Eq. (|2]). 
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